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Abstract 

Statistical mechanics of relative species abundance (RSA) patterns in biological networks is 
presented. The theory is based on multispecies replicator dynamics equivalent to the Lotka- 
Volterra equation, with diverse interspecies interactions. Various RSA patterns observed in nature 
are derived from a single parameter related to productivity or maturity of a community. The 
abundance distribution is formed like a widely observed left-skewed lognormal distribution. It 
is also found that the "canonical hypothesis" is supported in some parameter region where the 
typical RSA patterns are observed. As the model has a general form, the result can be applied 
to similar patterns in other complex biological networks, e.g. gene expression. 
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1 Macroscopic ecological patterns as eco- information 

The most significant feature of large-scale biological networks, such as food websmi2]Cmetabolic net- 
works in a cell[3J and protein networks^ is overwhelming diversity of components, e.g., species, 
chemical constituents and proteins, respectively, great complexity of network topology, and homeo- 
static stability of dynamics on the networks. From a theoretical viewpoint, it is a serious question 
how living organisms has evolved such a homeostasis because chaotic instability is inherited even in a 
simple nonlinear system. 

As an approach to such a problem, macroscopic patterns observed in various complex networks 
have been studied[H]. Such studies on scale-free networks have been elucidated the characteristics of 
topology of natural and artificial complex networks, and the evolutionary conditions which produces 
such a topology. Unifying approaches to biological and abiological networks have emphasized their 
similarity and difference. For example, it is pointed out that an infection of computer virus on the 
internet with scale-free topology is followed by a qualitatively different epidemic dynamics from the 
one of biological viruses in nature, and, therefore, the computer viruses are hardly eradicated 

On the other hand, in a large-scale complex biological networks such as ecosystems, not only 
a topology of the network links but also a thickness of each link, i.e. the strength of interactions, 
definitely affect population dynamics and resulting macroscopic patterns. In ecology, classical macro- 
scopic patterns observed and studied for a long time is RSA patterns, in other words, abundance 
distribution of species, which is one of the most accumulated informations obtained in ecology. 

Nevertheless, how to clarify the mechanisms underlying those RSA patterns has been one of the 
'unanswered questions in ecology in the last century [5]' even though the knowledge obtained from 
it would affect vast areas of nature conservation. Various models have been applied to ecosystem 
communities where species compete for niches on a trophic level E3 HH 1121 1131 1141 1151 1161 1171 
HH1 HU EOl EH 1221 EH ESI 123 CBl BHI, but these models have left the more complex systems 
a mystery. Such systems occur on multiple trophic levels and include various types of interspecies 
interactions, such as prey-predator relationships, mutualism, competition, and detritus food chains. 
Although RSA patterns are observed universally in nature, their essential parameters have not been 
fully clarified. In this paper, it is presented that RSA patterns are derived from a statistical mechanical 
theory [3D], based on a general evolutionary dynamics which is applied in vast area of fields. 
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We consider here the replicator equation 31 (RE), 



AT 
iV 



where iV is the number of species, and < Xi(t) < N denotes ith species' population. The functions 
fi(x) and f(x) denote fitness of species i and it's average, respectively. Interaction between ith species 
and j is specified by Jy . Note that total population is conserved at any time as ■ Xj = N and, that 
is, the trajectory of the dynamics is bounded in a simplex J^. ajj = iV. 

The RE appears in various fields 31 . In sociobiology, it is a game dynamical equation for the evo- 
lution of behavioral phenotypes; in macromolecular evolution, it is the basis of autocatalytic reaction 
networks (hypercycles); and in population genetics it is the continuous-time selection equation in the 
symmetric (J,j = Jj{) case. The symmetric RE also corresponds to a classical model of competitive 
community for resources |32j. The replicator dynamics, therefore, are often used as a model of complex 
systems in which many components changes their numbers through complex reaction, replication and 
reproduction of the components. 

Here we assume that {Jij) is a time-independent random symmetric (Jy = Jji) matrix whose 
elements have a normal distribution with mean m(> 0) and variance J 2 /N as 



N 

P(J - ) = ^ eXP 



— — ) (Jii — m ) 2 
2J 2 ' 



(i + j) 



Self-interactions are all set to a negative constant as Ja = —u(< 0). Note that the essential parameter 
is unique as p = (u + m)/ J because the transformation of the interaction Kij = (Jjj — m)/J does 
not change the trajectory of the dynamics (JIJ. Although ecologists do not generally believe in the 
randomness of interspecies interactions in nature, the discipline has been affected by the random 
interaction model |33j as a prototype of complex systems. 

Particularly in the context of ecology, the N species RE 31 is equivalent to the N — 1 species 
Lotka-Volterra (LV) equation 

dy 4 

= Vi I n - ( 



df 



JV-l 



That is, the abundance yi and the parameters in the corresponding LV are described by those in the 



3 



present RE model as, 



y i = x i /x M (* = 1,2,...,JV), (2) 

Ti = JiM — JMM = JiM + U, (3) 

bij = Jij — Jmj (4) 

where the 'resource' species M (i/m = 1) can be arbitrarily chosen from N species in the RE. The 
ecological interspecies interactions ^ j) have a normal distribution with mean and vari- 

ance 2J 2 /N from Eq. (gj, and they are no longer symmetric (6y ^ bji). The present model 
therefore describes an ecological community with complex prey-predator interactions ((bij,bji) — > 
(+, — ) or (— , +)), mutualism (+, +) and competition (— , — ). Moreover, a community can have a 'loop' 
(detritus) food chain — ► (+, -), {bjk,b k j) — > (+, -), {bki,b lk ) — > (+, -)). The intraspecific in- 
teraction &jj turns out to be related to the intrinsic growth rate r$ as &jj = Ju — J mi = — u — JiM = —i*,; 
and is therefore competitive (pa < 0) for producers (r^ > 0) or mutualistic (bu > 0) for consumers 

in < o). 

By Eq. the intrinsic growth rates also have a normal distribution with mean u+m and variance 
J 2 /N . The probability at which fj is positive-that is, that the i-th species is a producer-is therefore 
given by the error function, 

r°° Af 
Prob(r, > 0) = / —= exp (-i 2 ) 

Consequently, the parameter p can be termed as the 'productivity' of a community because the larger 
the p, the greater the number of producers. This can be also understood from the fact that p is 
connected to the average growth rate: 

■jr ^ n = {JjM + u)j = m + u=pJ. (5) 

i 

The parameter p is also connected to the maturity of an ecosystem because m increases in time in an 
evolutionary model |34|. 

Note that the growth rate (1 / Xi)dxi / dt in RE Q has no ecological meaning because it is defined 
by the average fitness / subtracted from the fitness fa. We therefore consider the equivalent LV 
when we discuss the model in the context of ecology as stated above. Why we do not consider LV 
with random asymmetric interactions from the beginning is that the system we consider here is not 
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general asymmetric LVs but a class of LVs which is corresponded to a symmetric RE whose symmetry 
is crucial to the present analysis. 

While random asymmetric interaction matrix ( ^ Jji and, Jy and Jji are independent each other) 
was assumed in the classical random population models j^S] 1331 I36| . here, the symmetric matrix 
(Jij — Jji) enables us to derive RSA patterns and left- skew ed^\ canonical\14l IT7| lognormal-like 
distribution, from a single parameter p|3U|. 



2 Mean field theory of random symmetric replicator dynamics 

The symmetry (Jy = Jji) makes the average fitness / = Ylj k JjkXjXk (the second term of the 
r.h.s. of Eq. l|T[l) a Lyapunov function (Appendix A) |31| . which is a nondecreasing function of time in 
dynamics l(T|l. Therefore, every initial state converges to a local maximum of / as t — > oo. Interpreting 
"K = — as an energy function, macroscopic (thermodynamic) functions of the system is derived 
from free energy 

f = — lim lim - „ T } J (6) 

J /J^ooA^oo N(3 

at such a maximum by using the technique of statistical mechanics of random systems [371 1381 1391 141)1 
EH Eg. The bracket 

f N oo \ 

" /«) (7) 



denotes the random averaae$7\ over random interactions, by which a typical behavior of the system 
can be analyzed. The normalization factor Z denotes a partition function with the condition (Y). Xi = 
N) and is represented as 



H dXi 5(N-Y, Xi )e-W=Tre-^, (8) 



where ensemble average is represented by the trace s\. 

Why we execute the random average is that free energy is "self-averaging" , that is, it is represented 
by an average over random interactions, not by a detail of each sample of interactions, which is justified 
in the thermodynamic limit, where the number of species are very large in the context of ecology. The 
inequality i < j in JjJ denotes the product of the values of? and j satisfying !<•••<£< j <-••<. N, 
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If we define Hamiltonian as 



J£ = -- ^ JjjXjXj + hy^^6(y- Xj), (9) 



2 

hi 



where 9(z)[= l(z > 0);0(z < 0)] is the step function, we will be able derive cumulative distribution 
function C (y) (proportion of the number of species which has abundance less than and equal to y) as 



C(y) = lim lim lim / -Tr fiA X J c -w\ = _ lim lim lim 



1 d(\nZ)j 

p'^So N^oo h"b \Zi Xi ^ N / j P^oa N^oo N (3 dh 

= } im %> ( 10 ) 

h^o ah 

where 8(x)[— l(x > 0);0(a; < 0)] is the step function. Information of equilibrium of RE QJ can be 
obtained by setting h = 0. By the identical equation 

(\nZ)j = lim ^ J - 1 (11) 

n-vO n 

the random average of the logarithm of the partition function, which is hard to execute analytically, 
can be transformed to the random average of an equivalent rt-replicated partition function 



<^ = ( n Sf~^ a > = ( n $} \-p - 2 e +/i E^- o 

\a=l / J \a=l V ij i 

which is more tractable. Each !K a denotes a replica Hamiltonian where an variable X4 is replaced 
by 1" with a replica index a — 1,2, ... ,n in Eq. (j5J). The analysis using above transformation of 
the Hamiltonian is called as the replica method 37 which enables us to analyze typical behavior of 
free energy with random interactions. The replica method was originally invented for analysis of 
the spin glass, magnetic alloy. Recently it has been successfully applied to various models 43| with 
time-invariant random interactions other than physical systems. By exchanging the order of the trace 
and the random average, we can precedently execute the random average. As the integrals (J7J are 
Gaussian integral of N(N — l)/2 variables Jy, the replicated partition function becomes 



(Z n )j = $[ } exp/3 



w E E x ^ + E E -lEE^-^EE^-o 

i<Cj \ a / a i<cj a % a % 



G 



As the first term in [• • • ] above can be rewritten as following, 



e -It 

a<b \ i I i 



E«) 2 



i 

+ 2 



-i 2 



E«) 2 



we can derive 



a<6 



" 1 4iV ^ 



E«) 2 



AN 2.^ 



-I 2 



E«) 



+ m££<^££(<) 2 -/ l ££%-<) 



a i<j 



if 2 



^3 



As we rewrite like following, 



a (. \ i / i 



K4. 



and apply the Hubbard-Stratonovich transformation 

I f°° fx 2 , \ 

-== / dx cxp h V 2Xax ) 

V 2 j 



^ _ 



where A > 0, the quadratic terms of Q2 i • • • ) in K\,Ki, K 4 can be transformed to linear terms. By 
this, we can execute the trace .Tj and obtain 





= n° xp 




a<b 


exp(/3X 2 ) 


= n° x p 

a 


exp(/3K 4 ) 


= n° x p 



(pj) 2 

2N 

m 2 

AN 



exp(,3A' 4 ) = n «p fe = n ^ /" * ^ |-f + V^E < ( " 
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If we transform the variables as y a b = Pv NY a b, s a = /?V NS a , t a = f3vNT a , we obtain 



(Z n ) — Tr 
V* h {x^} 



n 

.a<b 



dY ai 
L 



n 



dSg 

L 



n 



00 dTa 



x exp/3 



f^E^+^E^E^ 



V 2 

(f 



a<b 



E* 2 +(^)E^E«) 



(^)E^ 2 + v^E^E< 



a i 



ti + Til 



where L = tJ2tt/[3 2 N. Let us write g a b for the terms with J2i m [''']■ The delta function in the trace 
(JHJ can be represented as a Fourier transformation 



poo pioo 1 [ 



The terms including t/ Q {, and the terms of J^- above can be represented by a product of independent 
terms, and therefore can be written by the iVth power of a term in which index i is omitted as 



Tj } exp (^g ab - E r - E = (II J o °° dx ^ ex P U* - E r ^ a ) = CX P HA N = cxp N ln[A] 



The term g' ab denotes g a b without index i. Now we come to sum up the n- replicated partition function 
averaged over samples as 



n 

a<b 



dY ab 



n r £ 



nf ^ 



rr / dr. 



a JVG{F,S,T,r} 



(12) 



where 



G{y,5 > T,r} = -^- ^y^ + E^ + E^ 2 + + * 

\a<6 a a / a 



Tr 



exp Sab 



and = (J\ a dx a ) . By the saddle point method, the integral l|12fl can be replaced by the 
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integrand exp(iVG), and by substituting this for @) and Ijllfl . the free energy can be represented as 



1 



= — lim lim . 

Af^oo n^O \_ (3Nn 

= — lim lim 



Af^oo n^o Y f3Nn 
1 



= — lim lim 



iV->oon->o \f3Nn 

lim — Min/„ 

n— >o n 



{exp(-/37VMin/ n ) - 1} 

f rat v Min /" i , 
exp —pNn km | — 1 

\ n— >o n 



1 _ ^ i im , | 

n— >o n 



by the minimum of f n . If we transform the variables as q a b = Y a i,/J, cr a = S a \/2/ J, r a = T a ^J ' (3 / ! Ni 
the free energy becomes 



a<b 



/3iVm 



(13) 



where 



a<b 

u + m 



-f3 J 2 ]T ^ V ~ ^ >Z(x a ) 2 Va +NmJ2 * a Ta + ~ £ 

a a 

2 -}~2(x a r + h}~2e(y-x a ). 



Here we assume the Replica Symmetry (RS) as 

q = Qab, (r = <r a , t = T a , r = r n for Va, 6. 

By the discussion on the stability of saddle point solution for the estimation of the free energy, RS is 
justified at least for p > \/2 [381 I39| . By the RS order parameters, we can write as 



(/3J) 2 f _„ \ a ( I3J 2 _ (3.P _ u + m 



2 



5>° ^-^^ £(*•) 



(r + iV/3mr) - /%^] % ~ ^ a )- 



Bi 
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Here we can execute the trace in HKift and we obtain 



= II /°° d ^ a ex p I B i(E x<1 ) 2 + ^ I> a ) 2 - B * E - B A 

a J° [a a a J 

/•oo i ^oo ( i J 

= TT/ dx a -= dze- z2 / 2 cxp<^ x /2B' i y2x a z + B 2 y2(x a ) 2 - B 3 J2x a - B 4 \ 

a J V27T J-oo ^ la a a J 

y da; exp jzx^Bi + B 2 £ 2 - B 3 ie - /3h6(y - x)| 



dpi (2) 



At the final equality, we have replaced the rt-fold multiple integral of x a by a integral over a variable 
a; without the replica index a because each integral of x a is independent each other. Moreover, as we 
expect to take the limit n — > 0, using a n ~ 1 + nlnaAln(l + n) ~ n, we obtain 



In 



Tr p-^^^- 



In y dpi (2) 1 + nlnJ dxexp |s 2 .t 2 + (zy/2Bi - B 3 )a; - /3/i0(y - x) j 



In 



/> />oo /> />oo 

1 + n dpi (z) In / dxexp{---} ~n / dpi (z) In / dxexp{---} 



We, then, finally obtain the free energy density as 



/ = lim-Mm/„= limMini^(n-l)q 2 + ^a 2 - — 



r 1 

~ 



dpi (z) In / dare^^ 



Mini - &jLf + ^ 2 - ^r 2 -I * / dpi(z) In / dxe^ 
4 4 2 pp. 



(14) 



where 



<7(x) ee - j ^-g - ^-a + j x 2 + | z _ ( r /0 + JV"mr)| x - W(y - a?). 



Condition which minimizes /, which is the term in {} in the right side of Eq. (|14f) gives mean field 
equations. As the term with h is virtual for cumulative distribution function, we let it to be zero 
hereafter. By the calculations, df/dq = 0, df /da = 0, df/dr = 0, df/dr = 0, we obtain the 
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mean field equations, 



h^L H x ~m )— (15) 



a = J dpi(z) J dxx 



c /3ff(a) 
2 f 

Z 9 . 



(16) 



r r°° e Pg(x) 

J d Pl (z) J dxx^— = l ( IT) 



z 9 

where Z g = / °° dx exp(/3g(x)). By substituting (fTTTf) for l(TK|) . we obtain 

- g) = y d Pl (z) z ^ Zg PV ^ V " (18) 

To estimate the integral over x in the right side, let us rewrite like following 

g(x) = -\{{u + m)- J 2 f3(a-q))}x 2 + {J^/qz - {r / j3 + N mr)} x 
2 ^ ^ v > 

The condition E\ > should be satisfied for the convergence of the integral of a; in (|18fl . When 
E 2 /(2E 1 ) > 0, that is, z > (r//3 + Nm,T)/(J^/q) = -A in ifT^. the integral of x in Eq. (JTHJ) can be 
replaced by the integrand at the apex x\ = E2/2E1 in the limit (3 — > 00. We, therefore, obtain 

2E\ 



As 



Z g = / dxe^ = e ^/ iE \ 
Jo 



the contribution to Ijl8(l becomes E2/2E1. Similarly, when Eij2E\ < OAthat is, z < —A, the integral 
of x in Eq. Ijl8(l can be replaced by the integrand at X2 — 0, which is zero and has no contribution. 
After the similar calculation of the integral in Eq. Q17JI and some transformation of the expressions, 
we obtain the following mean field equations as 

/oo 
dpi(z)(z + A), 
-A 

\2 ' 



{p-vY = / d Pl (z)(z + AY 

J-A 

A = ^/q(p-2v), 
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where v = Jf3(a — q). As stated previously, the above result is essentially equivalent to the case 
m = 0, J = 1 |2H|. As the mean field equations can be solved analytically only for some special values 
of p (e.g., p = V%), we solve them numerically for a general range of p and obtain the order parameters 
9,u as a function of p, which are depicted in Figs, ^a) and 1(b). Macroscopic functions, such as 
diversity and abundance distributions, are simultaneously obtained by substituting q and v for them. 



3 Diversity and abundance distributions 

Among macroscopic functions calculated in the present framework, the most significant for a theory 
of RSA patterns is the cumulative distribution function of abundance (|10|l which is derived from the 
free energy / as 

C p (y) = \im^- = limj d Pl (z) J dx 9(y - x)eP 3 ^/Z g 

-A 



— A 



dp 1 (z)9(y-E 2 /(2E 1 )) + / d Pl (z)8(y), 



where the saddle point method was used in the last line like in the integral over x in HHJ. Substituting 
Ei and E 2 for the above and rewriting the population by x from ?/, the resulting cumulative distribution 
function is represented as 



C p (x) = C p (0)8(x) + / d Pl (z 



3G 



^q{z + A) 
x — - 



A V P~ v 

The quantity C p (0) = J_ A dpi(z) gives cumulative distribution function of species with zero popula- 
tion, that is, a proportion of extinct species. 

The function a p (0) = 1 — C p (0) = v(p — v) and a p (l) = 1 — C p (l) of p can be termed 'diversity', 
i.e., the proportion of nonextinct species and that of the species with abundance larger than unity, 
respectively, as depicted in Fig. [5J This demonstrates a typical positive correlation between produc- 
tivity and diversity Numerical results for a p (l) are also depicted in Fig.|2]for comparison. We 
see good agreement between the analytical and the numerical results for p > 1, while some deviations 
appear for small values of p. This small- value deviation is attributable to the occurrence of replica 
symmetry breaking (RSB) 37H3HI for P < which yields a number of metastable states of Eq. (|14|) . 
and the replicator dynamics essentially converges to not only a ground state of (|14fl but also to 
the metastable states. Since the energy "K and the diversity are both nonincreasing functions of time 
in dynamics Q), the mean- field results here give a lower minimum of diversity. Interestingly, the 
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metastable states enhance the diversity. The analysis of RSB is expected to improve the quantitative 
agreement 

In Fig[3 we also see a power law of the diversity S = Na p (Q) oc p v (p < 1; 77 ~ 2.3). This can be 
related to the species-area relationships S cx A x (A = const.) if p is a power function of area A, 

that is, larger the area, more producers are observed than consumers, which is one of the predictions 
in the present study. 

The function a p {x) = 1 — C p (x) is the survival function, the proportion of species whose abundance 
is larger than x. The survival function has been often used in the medical statistics. Note that a p {x) 
is also represented as a function of species rank n: 



a p (x) = ^ for i6[i (n+1) ,i W ), 

[n = 1,2, ... ,S < N) if the species abundance is ranked in descending order, as in x^ > x^ > 
■ ■ • > x^ > •• • > x^ > 0. As the function a p (x) is a nonincreasing monotonic function, the species 
abundance relation, i.e., the abundance x^ as a function of a rank n, is given by the inverse function 
of a p {x) as x^ — x p (n/N) — ap~ X ^[x), depicted in Fig. Inland Fig. 0]for some values of p. 

We observe two typical RSA patterns in different regions and with different species composi- 
tions ^5]: one is a straight line like the geometric series |^1 for a small value of p, and the other consists 
of sigmoid curves on a logarithmic vertical axis for some range of p. This latter RSA pattern denotes 
a lognormal-like abundance distribution. Remarkably, the transition of the RSA patterns from low 
p to high is identical to the observed transition from low- to high-productivity areas; that is, from 
a species-poor area such as an alpine or polar region to a species-rich tropical rain forest J21 • The 
transition also corresponds to the secular variation of RSA patterns observed in abandoned cultivated 
land |16| . This supports the contention that p (orm) is a maturity parameter, as is suggested by an 
evolutionary model |34j . 

The abundance distribution is also derived from the cumulative distribution function C p (x) by its 
derivative as 

" exp /_ bz*L L _ ifiLiMVl + Cp(0)J(l) , (20) 



y/2irq I 2q \ p — v 

and it is depicted in Fig. for some values of p. The first term is a normal distribution but not 
a lognormal distribution. Nevertheless, the curves in Fig. 01 demonstrate a typical sigmoid pattern 
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on a logarithmic vertical axis. This pattern indicates the coexistence of very abundant species with 
rare ones. This multiscale of abundance is intuitively understood by a divergent behavior of the 
variance a 2 = q/(j> — v) 2 of F p (x) for small p because q — > oo and v — > for p — > 0. Moreover, 
the mode of F p (x) per 'natural' octave Jl] ln(x) is always a positive value (as shown in Fig. at 
x* = j (A + VA 2 + 4) > 0, which denotes a unimodal distribution. Indeed, the mode diverges as 

„ _a_ = 1 > 

~* |A| ~ (p-v)\p-2v\ 

for p — > 0. As a result, the abundance distribution is a normal distribution truncated at x = and 
given in the positive abundance range x > and it has a large variance <r 2 — > oo and a negatively 
divergent mean /z = — oo satisfying £ = ^ — * for p — * 0. This indicates that, for small 

p, Fp(x) becomes a tail of broad distribution but it still has a peak at a positive abundance x* when 
plotted on the log-scale horizontal axis. This is why the abundance distribution per octave looks like 
a left-skewed lognormal distribution ^5] in Fig. [S] 



4 Canonical hypothesis 

According to the canonical hypothesis |14l I17 j . the quantity 7 = log(a;jv)/log(a; max ) takes a value 
near unity in various real communities, where x max is the abundance of the most abundant species 
and Xn gives the position of the mode of individual curve P p (x) = xF p (x). Using the abundance 
distribution (|20ll . we derive an analytical expression for the above functions to check the validity of 
the canonical hypothesis. We first evaluate an expected value of the most abundant species x max . 
From the definition of x max , that is NF p (x ma x) — 1, and the conservation of the total abundance 
J °° xF p (x)dx = N, which is equivalent to X{ = N, we obtain 

q(p -v) + a^2\n(^ ^ Aa ^ 0) ) 

Xmax 

p — V 

On the other hand, the mode of the individual curve P p (x) per octave is given by xn — § ( A + VA 2 + 8) , 
and finally, the parameter 7 is evaluated by substituting the values of the order parameters q and v for 
each value of p. In the present model, 7 is a monotonically increasing function of p and 0.96 < 7 < 1.04 
for 0.1 < p < 0.6, denoting that the canonical hypothesis is supported in the range of p giving the 
typical RSA patterns in Fig. El Although the canonical hypothesis was demonstrated to be merely a 
mathematical consequence of lognormal distribution |17| rather than anything biological, it is notewor- 
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thy that the lognormal-like abundance distribution with 7 ~ 1 derives from basic ecological dynamics. 
This still suggests a biological foundation for the hypothesis in a large complex ecosystem, in the same 
way that a biological foundation was indicated for the theory of a local competitive community |18] . 

5 Topology of interactions 

The present theory seeks to capture the influence of productivity on the RSA patterns under the 
assumption that all species interact randomly; nevertheless, this assumption itself is never justified 
because it ignores a biological correlation between interactions produced by evolution. However, note 
that the randomness is assumed only for an initial state with N species in Eq. JJJ. Actually, the 
simulation reveals the resulting interactions of nonextinct species to be nonrandom. 

In Fig. E| interspecies interactions (bij) between the non-extinct species of the corresponding LV 
equations Q are depicted, which is obtained by numerical simulations of RE Q and the transfor- 
mation © and (J3J. The numerical integration of RE was executed for initial diversity N — 2048, a 
randomly generated interspecies interactions (Jij), a random initial population (xi), and for (a) p = 0.1 
and (b) 0.2 by the fourth-order Runge-Kutta method. At the equilibrium of Q for the parameters, the 
number of non-extinct RE species Na p (x = 1) was (a) 11 and (b) 28. In each figure, the non-extinct 
LV species (yi > 0) without the resource species (%)m = 1) are depicted by a blue disk which is arranged 
clockwise in descending order of the intrinsic growth rate r» as 1 > n > T2 > ■ ■ • > r, > • • • > rt > 
(hence every non-extinct LV species is a producer) where the number of non-extinct LV species is (a) 
L = 10 and (b) 27. The diameter of the disk is in proportion to | log(ri)|/| log(rj)|. Each type of 
interaction is represented by its color: green links denote a mutualistic interaction (6^, bji) = (+, +), 
yellow, a competitive (— , — ) and blue, an exploitation of more productive i on less productive j(> i) 
(+, — ) for i < j. No exploitation of less productive j on more productive i{< j) (— , +) is observed (it 
were drawn by a red link). The thickness of each link is proportional to the larger value of \bij\ and 

IM- 

It should be noted that not only each sample (a) and (b) but also every sample calculated for Fig. [21 
evolved to only flora, Vi Ti > 0. In every sample, only green, yellow and blue links are observed but no 
red; thus the stable community after extinction dynamics obtains a hierarchical structure. Moreover, it 
is observed that more productive species with larger tend to have competitive (yellow) relationships 
each other and less productive ones have mutualistic ones (green), the quantitative estimations of 
which is now in progress. It is suggested that this emergent hierarchy is connected to the stability of 
a large-scale plant community with complex interspecies interactions of not only competition but also 
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mutualism and exploitation. 



6 Discussions 

In the present model, all species coexist only in the limit p — > oo, that is, in the trivial cases in 
which interspecies interactions are negligible (J <C u ) or homogeneous (J — > 0), thereby giving 
otoa{x) = 0(1 — x), x' n ) = Xoo(n/iV) = 1 for all n and F oa (x) — 5(x — 1). Such a region of too 
large productivity, therefore, never corresponds with a real community even if all species coexist. The 
symmetric interactions (Jy = Jji) considered in the present study really emerged in a evolutionary 
community-assembly mo del [51]. In the simulation, though the introduced mutants had asymmetric 
interactions with existing species, the system evolved to have symmetric interactions. The system 
also showed a typical RSA pattern and therefore it can be said that the present study is an analytical 
treatment for it. The system moreover was resistant to exotic species, which is manifested in the 
functional form of v in Fig. ^b). As the order parameter v corresponds to susceptibility to external 
noise in the context of statistical physics, Fig. ^b) suggests that an ecosystem with medium produc- 
tivity (p ~ V2) is more sensitive to external disturbance than ones with lower or higher productivity. 
In other words, in the range p < \/2 where the typical RSA patterns are observed, the lower p, the 
RSA patterns are more robust to external noise such as environmental change, which is one of the 
predictions of the present study and can be verified by field studies. 

Simplicity of the present model with only one parameter conduces to some predictions to be verified 
by experimental researches. They are summarized as follows. 

(1) A RSA pattern of a community with not only competition but also mutualism and exploitation 

shows itself like a tail of broad normal distribution. Although such a truncated normal distri- 
bution with negatively large mean and large variance has not been examined to fit field data, 
there is still plenty of room to consider alternative distributions for a community with complex 
interactions, to which the other models of competition, such as niche apportion models [TT1 12()j 
or the neutral model |23> is n °t applied. 

(2) Diversity (number of non-extinct species) S is a power function of the productivity p which is 

proportional to the average growth rate in Eq. (JSJ. As suggested in the first paragraph of 
this section, the present theory is not valid for a too large value of the productivity, and the 
power law here may be applied to the left half of the hamp-shaped relationships with a peak at 
intermediate productivity levels, which has been reported most widely in field and experimental 
researches 013 El • 
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(3) Productivity (the average growth rate) p is a power function of area A. If it is, the present 

theory also predicts the species-area relationships. The dependence of p on A appears against 
the intuition that the average growth rate per species is constant even if we enlarge an area 
of observation. It should be, however, noted that such a constant p is justified only if the 
distribution F p (x) is invariable under changes of area A. Although the present analysis assumes 
infinite total population N, thereby infinite area A, the finite-area effect on p and F p (x) will be 
an important subject to be addressed in future studies. 

(4) The transitions of the RSA patterns from species-poor to species-rich community or from im- 

mature to mature community attribute to the productivity or the average growth rate p. Such 
relationships between the various RSA patterns and ecological parameters will be one of the 
focal points of the next generation of community ecology though some classical models have no 
parameter and often gives no explanation on variations of RSA patterns. 

(5) The canonical hypothesis is supported. Moreover, the value 7 is increasing function of p. Com- 

pared to the distribution F p {x) itself, the statistical evaluation of which is often controversial 23 , 
quantities like 7 seem to be more tractable in quantitative study of field data and there still 
is a room for consideration of such macroscopic quantities which may characterize a large and 
complex community. 

(6) A stable and complex community has a hierarchical structure in which more productive species 

exploits less ones, more productive species compete each other, and less productive species have 
mutualistic relationships among themselves. Exploring such a hierarchy and the bias of the 
competition and the mutualism will be one of the clues to clarify the unsolved problems on the 
complexity and the stability in community ecology. 

Verification of the above predictions are in progress through collaborations with field ecologists. 

In summary, it has been demonstrated that empirically supported patterns are derived from a single 
parameter of general population dynamics. This not only suggests the importance of globally coupled 
biological interactions in a large assemblage but also provides a unified viewpoint on mechanisms 
of similar patterns observed in other biological networks with complex interactions; for example, a 
lognormal abundance distribution of a protein in cells [48 ,49,50, which is revealed by gene expression 
networks. 
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A Average fitness as a Lyapunov function 

Since the interaction matrix J is symmetric (Jjj = Jji), the time derivative is written as 



^-(l)<-* >-(£)■'■-' (£)-*(£'- 



It is, therefore, found that 



I d/ft) 
2 di 



dx ? ; 



= ^^(Ja;) 2 - 



Xi {(Jx)i — x ■ Jx] 2 > 0, 



and the average fitness / is non-decreasing function of time, a Lyapunov function. 
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Figure Captions 

Figure ^ Order parameters (a) q and (b) v as a function of the productivity parameter p. 

Figure |21 Diversity a p (x — 0, 1) as a function of p of log-log scales for x — 0(red) and x = l(green). 
Black circles show numerical solutions of a p (l) averaged over randomly generated 50 samples 
of (J i:j ) for Eq. Q with N = 2048 and p = 0.1, 0.2, 0.3, 0.4, 0.5, a/2/2, 1, V%, 2, 3. Numerical 
integration of RE was executed by the fourth-order Runge-Kutta method. Error bars indicate 
the maximum and minimum values found in the samples. 

Figure [3] Rank-abundance relations as a function of productivity p on normal-normal scales. 

Figure [3] Rank-abundance relations as a function of productivity p on normal-log scales. 

Figure |5] Abundance distribution per 'natural' octave ln(a;). Functions F p (x)x(d(ln(x)) — F p (x)dx) 
for some values of p are depicted, whereas Preston originally defined octaves as logarithms to 
base 2 [Hj. 

Figure |6] Interspecies interactions (fry) between the non-extinct species of the corresponding LV 
equations. 



23 




0.5 1 1.5 2 2.5 3 



P 

Figure 1: (a) TOKITA 



24 




0.5 1 1.5 2 2.5 3 

P 



Figure 1: (b) TOKITA 
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Figure 2: TOKITA 
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Figure 3: TOKITA 



27 




Figure 4: TOKITA 
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Figure 5: TOKITA 
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